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ABSTRACT 


Two  dynamic  methods  are  compared  for  integrating  the 
Reynolds  equation  as  applied  to  determining  the  axial  dynamics 
of  a  spool  bearing.  It  is  shown  that  very  sensitive  phase 
shifts  in  the  numerical  schemes  can  badly  falsify  the  dynamics 
and  thus  care  is  needed  in  interpreting  results  using  classical 
methods.  Second  order  temporal  accuracy  can  alleviate  these 
problems  and  one  such  scheme  is  presented  with  excellent  phase 
properties  that  is  just  as  computationally  efficient  per  time 
step  as  any  of  the  first  order  accurate  schemes. 
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Introduction 


A  large  number  of  papers  have  appeared  in  the  gas  bearing  literature 
solving  the  time  dependent  compressible  Reynolds  equation  by  direct  numerical 
integration.  A  typical  example  of  modern  practice  is  given  in  [1].  In  some 
of  our  recent  work  a  notable  example  of  the  sensitivity  of  time  transient 
dynamics  to  the  accuracy  of  the  integration  scheme  appeared.  The  most  im¬ 
portant  result  was  the  failure  of  the  classical  implicit  scheme  to  adequately 
track  the  system  dynamics  except  when  the  time  step  was  very  small.  It  will 
be  shown  that  this  defect  is  due  to  phase  shifts  in  the  numerical  solution 
and  that  these  problems  can  be  adequately  cured  by  using  an  extrapolated 
Crank-Nicolson  scheme. 

Such  phase  shifts  can  have  catastrophic  effects  in  the  numerical  simu¬ 
lation;  however,  they  are  dependent  on  the  particular  integration  scheme 
used.  The  fact  that  a  scheme  is  inherently  0(aT)  in  itself  is  a  priori 
insufficient  to  categroize  the  scheme  as  having  poor  phase  properties.  For 
example,  in  the  classical  work  [2]  the  infinite  length  journal  bearing  was 
solved  using  PH  as  the  dependent  variable  (a  natural  choice).  In  some  of 
our  own  unpublished  work,  the  calculations  of  [2]  are  corroborated;  further, 
however,  both  the  schemes  discussed  here  are  used  (solving  for  PH),  and 
phase  shifts  are  highly  mollified,  although  the  extrapolated  Crank-Nicolson 
schemes  proves  itself  the  superior.  In  diverse  applications  such  as  [1], 
careful  a  posteriori  testing  is  necessary  to  validate  the  system  behavior 
extracted  from  the  numerics. 
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The  problem  of  interest  was  to  solve  for  the  axial  dynamics  of  a 
spool  bearing  with  spiral  groove  thrust  plates  and  herringbone  groove 
journal  (See  [3],  Figures  1  and  2).  With  no  radial  eccentricity,  the 
geometry  is  axisymmetric,  and  thus  spatially  the  problem  is  one  dimensional. 
Appendix  1  depicts  the  equations  of  interest  according  to  the  narrow  groove 
theory  and  Appendix  2  derives  the  equations  of  motion  for  the  axial 
dynamics. 


Implicit  Scheme  and  Results 

In  Appendix  3  the  basic  spatial  discretization  for  the  divergence 
form  of  Reynolds  equation  is  derived.  The  implicit  scheme  approximates 
the  highest  spatial  derivatives  at  the  advanced  time  level.  All  other 
terms  are  evaluated  at  the  current  time  level.  Thus  (33)  takes  the  form 


n  .  _rrl-l  n 
aJ  +  .5AJ  ^  +.5 


-  a 
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-  1 
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where  superscript  n  refers  to  current  time  level,  and  n+i  the  advanced  time 

level.  As  discussed  in  Appendix  3,  the  last  term  on  the  right  only  appears 

at  the  chamfer  and  if  there  is  no  chamfer  C.  =  0. 

n 


a) 


The  equations  of  motion  (28)  given  in  Appendix  2  are  integrated  by  a 
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standard  Euler  scheme  accurate  to  0(aT),  which  is  consistent  with  (1); 
thus 


(if1  '  if /AT  =  l* 

(ef 1  -  ef  /AT  =  (-a11  +  Fn/m)/c  v2 
i  i  e  o 


(2) 


(1)  and  (2)  can  be  computed  independently  in  parallel  to  obtain  the  new 
pressure  Pn+*  and  the  new  eccentricity  ef**  .  The  new  gap  is  then 


i  i  -n+1 
i+  ct 


(3) 


and  the  gas  bearing  force  at  the  advanced  time  is  obtained  by  integrating 
Pn+1. 


There  are  many  allowable  variants  of  (1);  e.g.,  evaluating  b  at  n+1 
instead  of  n,  but  these  are  still  0(aT).  The  choice  of  (1)  as  the  implicit 
variant  is  dictated  by  the  ellipticity  of  (20)  which  implies  the  tridiagonal 
matrix  to  be  inverted  to  solve  (1)  is  diagonally  dominant,  „nd  therefore 
Gauisian  elimination  without  pivoting  is  stable,  so  (1)  can  be  very  efficiently 
solved  (see  the  discussion  in  Appendix  A  of  [4]). 

The  spool  bearing  with  the  geometry  given  in  Appendix  4  was  tested 
for  stability  using  the  implicit  scheme  (1)  and  (2).  The  gap  was  given  an 
initial  perturbation  of  approximately  1.13  in. /sec.  Time  was  nondimensional - 
ized  by  letting  v  =  400  Hertz,  and  the  natural  frequency  of  this  particular 
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bearing  is  about  6400  Hertz.  In  Figure  1  the  nondimensional  eccentricity 
versus  nondimensional  time  is  shown  for  three  different  time  steps.  The 
abscissa  is  shown  in  a  new  nondimensional  time  T'=  T/(2-rt/288) ,  the  ordinate 
is  the  negative  of  the  nondimensional  eccentricity,  and  N*  corresponds 
very  closely  to  the  number  of  time  steps  per  cycle  of  natural  frequency. 

The  clear  picture  that  emerges  from  Figure  1  is  the  large  variation 
of  bearing  response  as  a  function  of  time  step,  and  even  as  Figure  2  shows 
with  N*  =  144,  the  bearing  appears  to  be  neutrally  stable.  It  can  be 
shown  that  this  bearing  should  be  stable,  so  that  the  case  N*  =144  still 
depicts  the  dynamics  incorrectly.  Such  a  posteriori  examination  of  the 
data  shows  the  time  step  must  be  much  smaller  yet  to  obtain  correct  results. 

A  slightly  deeper  investigation  explains  the  apparent  instability 
shown  in  Figure  1.  Let  us  consider  the  bearing  response  to  a  forced-gap 
oscillation.  The  equations  of  motion  are  uncoupled,  and  the  bearing  gap 
is  forced  by 


H  =  1  +  c* cos(T)  (4) 

where  e  generally  is  very  small  (in  this  case  T  =  .01).  The  bearing  force 
is  decomposed  into  a  Fourier  series  over  one  cycle 


F(T)  =  a 

O 


cos(iT)  + 


(5) 


Eventually  the  force  F  will  entrain  itself  into  a  periodic  function,  so 
the  coefficients  {a.}  and  {b.}  will  not  change  from  cycle  to  cycle. 

In  fact,  with  T  small,  only  a-j  and  b-j  are  nonzero.  By  the  usual  conventions, 
dimensionally  the  in-phase  and  out-of-phase  components  of  the  bearing  stiff- 
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ness  are  U  =  -a1/«cQ  and  V  =  respectively.  Table  1  shows  the 


resul ts 

for  v=  400  Hertz, 

for  N*  taking  on 

the  same  values 

as  in  Figure  1 

The  two 

columns  under  0(aT)  show  the  in-phase  and  out-of-phase  components 

for  this 

implicit  scheme. 

1)  and  V  are  in 

units  of  Ib/in. 

t 

TABLE  1 

.  U  and  V  for 

v=  400 

OUT) 

OUT2) 

N* 

11 

V 

U 

V 

18 

389,767 

-134,424 

440,680 

-42,782 

36 

414,829 

-  89,517 

437,861 

-41,074 

72 

426,269 

-  65,556 

— 

— 

X 

437,709 

-  41,595 

436,921 

-40,504 

LINEAR 

438,521 

-  39,810 

438,521 

-39,810 

The  row  marked  LINEAR  shows  the  results  obtained  from  a  small  pertur¬ 
bation  analysis.  The  major  source  of  error  is  the  large  shifts  in  V  as 
a  function  of  N*.  It  is  interesting  to  note  that  both  U  and  V  are  very 
linear  functions  of  aT  (equivalently  1/N*).  The  row  marked  X  is  the  result 

of  extrapolating  the  N*  =  36  and  N*  =  72  solutions  to  infinity  to  yield  a 

2 

solution  of  accuracy  0(aT  }.  It  is  more  descriptive  to  recast  F  in  the 
form 


F  =  Kcos  (T-<t>)  (-7c  ) 

O 

2  2  1/2 

where  K  =  (u  +  V  )  '  and  $  =  arctan  ( -V/U ) .  These  results  are  shown  in 
Table  2  with  <j>  actually  given  in  degrees  and  not  radians. 
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TABLE  2. 

K  and  *  for  v 

=  400  Hertz 

0(AT) 

0(aT2) 

N* 

K 

K 

18 

412,296 

19.03 

442,752 

5.54 

36 

424,378 

12.18 

439,783 

5.36 

72 

431,280 

8.74 

— 

-- 

X 

439,681 

5.43 

438,794 

5.30 

LINEAR 

440,324 

5.18 

440,324 

5.19 

The  major  source  of  error  can  now  be  seen  to  be  in  the  phase  lag 
<j>,  and  just  a  few  degrees  shift  cause  large  errors  in  V.  The  implicit 
scheme  induces  excessive  phase  lag  in  the  numerical  solution.  In  order 
to  see  why  the  excessively  poor  bearing  stability  properties  of  Figure  1 
arise  consider  Table  3  which  is  similar  to  Table  1,  but  v  =  7000  Hertz. 


TABLE  3.  U  and  V  for  v=  7000  Hertz 


0(at)  0(at2) 


N* 

U 

V 

U 

V 

18 

468,700 

-56,042 

476,473 

30,766 

36 

472,580 

-12,076 

474,425 

31,519 

72 

473,508 

9,699 

— 

-- 

X 

474,436 

31,474 

473,742 

31,770 

LINEAR 

477,575 

34,314 

477,575 

34,314 

7. 
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The  same  comments  apply  here  as  apply  to  Table  1.  For  N*  =  18  and 
36  case  V  shows  a  lagging  phase  rather  than  the  leading  phase  of  the 
actual  solution.  Again  it  is  important  to  recognize  that  V  can  be  very 
accurately  described  by  V  =  VQ  +  V^/N*.  and  it  is  the  excessive  size  of 
V-|,  a  characteristic  of  this  implicit  scheme,  that  degrades  the  results. 
From  Table  1  and  Table  3  we  are  able  to  infer  that  this  scheme  causes  a 
large  shift  in  the  negative  direction  of  V. 


For  one  degree  of  freedom  systems  the  solution  of  the  characteristic 
polynomial  is  given  by 


mv  =  U(v) 


V(v)=  0. 


(7) 


If  U  and  V  are  obtained  numerically  as  a  function  of  N*,  then  U  =  U(v,N*) 
and  V  =  V(v,  N*).  In  the  above  example,  V(v,N*)  is  approximately  in¬ 
variant  with  N*  and  always  U(v,  N*)  >0.  V(v,N*)  shows  large  shifts  in 

the  negative  direction  as  N*  goes  to  0.  The  zero  crossing  occurs  at 
larger  values  of  v(V  hasmonotone  increasing  behavior  as  a  function  of  v 
for  this  bearing  away  from  v=  0).  Since  1)  is  approximately  invariant 
with  N*.  by  (7)  m,  the  critical  mass,must  get  small  with  N*.  Thus,  the 
bearing  has  a  small  apparent  critical  mass,  which  explains  the  instability 
apparent  in  Figure  1  (which  worsens  with  smaller  N*). 


Second  Order  Scheme  and  Comparisons 


It  will  now  be  shown  how  an  extrapolated  Crank-Nicolson  scheme  cures 
the  preceding  defects  with  no  essential  penalty  in  computing  time.  The 
efficiency  of  this  scheme  was  shown  in  [4].  Define 


p1*"" 
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(8) 


The  extrapolated  Crank-Nicolson  scheme  is 


-.5  AJ  L'5(P  +  p  )J  ~bj+.5  '  aj- .5  AJ-1  L*5(P  +  P  )J  +  bj.  < 


pn+1  n 

_rrf.5\  J  _J  ,  r/.c  >-5,,.  _tr+- . 5\ 


.I  .  wf.5.  .  _  rrf.S'V  J  J.  „  n+.5._  n+.5\ 

-•5Wj+  +asj-i8j-  )—zr-+  •5Wj+  +asj-icj-  ) 


+  a—  —7—  r.c  (*?’5  [aH^i^/aT  +  aHttf:5/aT]  +  c  ) 

tq  oi.+a^  j  b  \  J  L  r,J+  r,J-  J  n  AT  J 


_n+l  _n 
P.  -P, 


Again  the  last  term  on  the  right  appears  only  at  the  chamfer.  In  any 
coefficient  evaluated  at  n  +  .5,  P  is  obtained  from  (8).  Hr  *  for  these 
coefficients  is  determined  as  shown  below.  Note  that  (8)  cannot  be  used 
at  the  initial  value  of  time  or  whenever  there  is  a  discontinuity  in  the 
equations  of  motion  as  a  function  of  time.  Whenever  this  occurs  (1)  is 
used  with  aT  changed  to  ,5aT.  Since  (9)  is  a  linear  system  for  Pn+\ 
it  is  no  more  time  consuming  to  solve  than  (1). 
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The  equations  of  motion  cannot  be  integrated  in  the  same  fashion 

2 

as  in  (2).  A  modified  Euler  integration  is  used  to  maintain  an  0(aT  ) 
accuracy.  This  is  a  two  step  integration  with  the  first  step  given  by 
(2)  with  aT  replaced  by  .5aT. 


n+.5 
e 


e”)/ .5 AT  =  ej 


(10) 


(e^"*5  -  e”)/-5AT  =  (-a^+  Fn/m)/cQv2 

To  integrate  to  the  full  step  it  is  necessary  to  obtain  the  gas  bearing 

71+  5 

force  F  at  T  +  .5aT.  This  is  accomplished  by  integrating  P  '  given 
by  (8).  The  second  step  in  the  integration  is 


.-rri-l  -n.  -n+.5 

(«1  -  eim  =  e2 

(11) 

,-rrf-l  -n.,_  .  n+.5  ,  „n+.5 ,  v.  2 

(«2  -  e2)/AT  -  (-ae  +  F  /m)/coV 


which  is  centered  about  n  +  .5.  Now  analogous  to  (3), 


1  +  .5 


,-n  -rri-1. 

(*L  +  )• 


(12) 


The  same  calculation  shown  in  Figure  1  for  N*=18  is  depicted  in 

Figure  3.  If  N*  is  repeatedly  doubled  there  is  negligible  shift  in  the 

results.  Further  contrast  in  the  quality  can  be  seen  in  the  tables. 

2 

Under  the  heading  0(&T  )  are  listed  the  results  for  this  scheme.  The  row 

3 

marked  X  is  a  quadratic  extrapolation  yielding  an  0(aT  )  error.  The 
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N*  =  18  solution  is  better  than  any  of  the  implicit  cases  shown,  and 
from  Table  2  it  can  be  seen  there  is  negligible  phase  shift.  Table  3 
indicates  both  U  and  V  are  accurate  at  the  higher  frequencies,  and 
thus  from  Tables  1  and  3  there  is  no  shift  in  V  as  illustrated  for  the 
implicit  scheme.  This  explains  the  very  good  results  depcited  in 
Figure  3. 

It  is  interesting  to  note  in  the  tables  the  slight  systematic 
differences  between  the  extrapolated  solutions  and  the  small  perturbation 
analysis.  The  perturbation  results  were  obtained  from  another  program 
in  which  the  difference  approximations  were  derived  in  a  slightly 
different  fashion.  Consistent  use  of  the  same  difference  approximations 
would  yield  identical  results. 


Discussion  and  Conclusions 


The  proven  convenience  of  dynamic  simulation  programsto  validate 
bearing  designs  makes  these  programs  valuable  tools.  We  have  illustrated 
one  of  the  difficulties  that  require  careful  consideration  in  interpreting 
stability  results,  for  example,  in  complicated  dynamic  simulations.  Very 
fortunately  simple  a  posteriori  comparisons  can  elucidate  whether  numerical 
difficulties  exist  (although  not  necessarily  cure  them). 

It  was  shown  that  the  classical  implicit  scheme  with  an  0(aT)  truncation 
error,  because  of  the  size  of  that  truncation  error,  will  manifest  slight 
phase  shifts  in  the  numerical  solution  which  will  appear  as  errors 
particularly  in  the  out-of-phase  component  of  the  bearing  force.  This  re¬ 
sults  in  falsification  of  the  stability  behavior  of  the  bearing.  An  extra- 

2 

polated  Crank-Nicolson  scheme  with  0(aT  )  truncation  error  cures  these 
difficulties  and  has  excellent  phase  properties. 
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What  was  not  considered  was  the  determination  of  optimal  0(aT) 

schemes  which  minimize  particularly  phase  errors.  Thus,  the  possibility 

is  not  ruled  out  that  variants  of  an  implicit  scheme  will  have  much 

better  phase  properties,  however,  it  is  unlikely  that  such  schemes  can 

2 

match  the  excellent  accuracy  of  the  linear  0(aT  )  scheme  demonstrated 
here. 


Changing  bearing  parameters  is  not  expected  to  make  significant 
changes  in  the  results.  When  A  was  reduced  by  a  factor  of  five  by  changing 
pa,  the  same  very  poor  out-of-phase  characteristics  depicted  in  Table  1 

a 

resulted. 

In  general,  if  testing  of  a  classical  scheme  indicates  numerical 
difficulties  (albeit  numerical  stability  or  truncation  error  effects), 
then  a  simple  alteration  to  an  extrapolated  Crank-Nicolson  scheme  can 
cure  both  problems  without  essential  penalty  in  computation  time. 
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parisons  were  supported  by  the  Office  of  Naval  Research  under  contract 
N00014-70-C-01 35  administered  by  Mr.  Stanley  Doroff. 
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APPENDIX  1.  Equations  and  Boundary  Conditions 

The  purpose  of  this  appendix  is  to  display  the  particular  one  dimensional 
parabolic  equation  being  solved  along  with  the  appropriate  boundary  conditions. 
We  are  interested  in  the  pressure  profile  in  a  spool  bearing  with  no  radial 
eccentricity  where  the  thrust  plates  are  spiral  grooved  and  the  journal  bearing 
herringbone  grooved. 

It  has  been  shown  [5]  that  according  to  the  narrow  groove  theory  the  one 
dimensional  Reynolds  equation  in  generalized  coordinates  for  all  spiral  groove 
geometries  in  nondimens ional  form  is 

?v  a 

-  5S--*V^2  **■  "W  ■  0  <13> 

where  the  mass  flux  is 

Ys  -  -PH3  KXR  ||  +  A6K4R2  cos  0  P  (14) 


A  mf(  fn)  (u,l+u>2> 


u>2*U>i 
A6  A  gu  2'to  ^ 


or  (1  -  a)  sin  0  ^ 


of  f  +  1  -  or 


0(1  *  sin  20  (r3  -  D2  +  r3]/  [(1  -  <*)  r3  +  a] 


1]  /[  (1  -  aor  +  Of] 


T  -  vt 


S  ■  s/r 


H  *  h  /c 
r  r  o 


•tf>  v  tflWr-siifr’v  ••  ■ 
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The  function  R(S)  is  the  normalized  bearing  radius  at  S;  if  Tj  is  the 

radius  of  the  journal  bearing,  then  R(S)  «*r./r  for  all  S  located  along  the 

j  o 

journal. 

All  other  parameters  are  conventional  spiral  groove  quantities  and  are 
defined  in  the  nomenclature  and  in  [  5  ]  . 

Let  us  now  consider  the  boundary  conditions.  At  the  ambient  edges  of  the 
thrust  plates  P  «*  1.  A  chamfer  region  may  exist  between  the  thrust  plates  and 
the  journal  bearing.  Two  possible  boundary  conditions  are  allowed  at  the  chamfer  - 
it  may  be  vented  to  some  arbitrary  pressure  Pc  which  need  not  be  Pq  =  1,  or  if 
unvented  then  it  can  be  shown  [5  ,  p.  28]  that 


out 


_§iL 


2  3 

Tip  C 

ra  o 


(16) 


where  Vc*  chamfer  volume.  If  the  chamfer  cross  section  is  shaped  as  an 

isosceles  triangle  with  height  c^  and  base  2c^,  then  the  approximate  chamfer  volume 

is  2rTr.c2  however  a  better  approximation  is 
j  h 


V  2TTrJch(ch  +  ht) 


(17) 


where  ht=*  thrust  plate  gap  it  the  chamfer.  (1?)  includes  the  chamfer  volume 
contribution  from  the  thrust  plate  gap  which  is  only  important  when  axial  time 
dependent  effects  are  being  considered.  Using  (17),  (16)  may  now  be  rewritten 


-y 

in  XS 


out 


c 

-Ar 

c 


2  V 


tu  1  +  ou2 


R.C.  j" P  +  C,  ||  "} 

j  h  L  ST  h  ST  J 


(18) 


where  C.  ■  c,  /c  and  R.  *  r./r  .  If  C,  =  0,  then  (18)  is  still  valid  as  the  interface 
n  h  o  j  j  o  h 

boundary  condition. 

With  a  little  algebra  it  can  be  shown  that 


&H 


~  (RK  PH  )  «  R  (k  H  ll  +  P  —  ) 
8T  u)  r  \  ®  r  ST  ST  / 

so  the  form  of  (13)  of  interest  is 


as 


-  A 


+  u)2 


R  (k  H 

\  (C  r  ST  9T  / 


(19) 


0. 


(20) 


APPENDIX  2.  Equations  of  Motion 

In  this  appendix  we  obtain  the  equations  of  motion  for  the  axial 
dynamics  of  a  spool  bearing.  Consider  an  inertial  frame  located  at  the 
center  of  the  journal’s  initial  position.  Let  xfc  be  the  displacement  of  the 
stator  and  let  it  have  mass  M.  Let  x^  be  the  displacement  of  the  rotor  and 
presume  it  has  mass  m.  Let  F  be  the  gas  film  force  and  f  the  external  force, 
then 

Md2xfc/dt2  **  f  -  F  (21 ) 

2  2 

md  x^/dt  ■  F  (22) 

combining  (21)  and  (22)  gives 

d2(xt  -  x^J/dt2  *  f  /M  -  F  (m  +  M)/mM  (23) 

Under  the  assumption  M  »nl,  m/ (m  +  M)  -m/M,  so  (23)  becomes 

md2  (x  -  x  )/dt2  *  mf/M  -  F  (24) 

t  U) 

us ing  the  convention  that  cq  is  the  nominal  thrust  plate  clearance,  it  is 
convenient  to  redefine 

X  -  x  -  c  -  e  (25) 

where  e  is  defined  as  the  displacement  of  the  rotor  relative  to  the  stator 
(which  is  the  dimensional  axial  eccentricity)  with  the  convention  that  e  is 
positive  in  the  direction  of  increasing  S  which  generally  would  be  the  direction 
of  the  positive  spin  axis.  By  definition  the  positive  direction  of  S  along  the 
journal  points  to  the  positive  thrust  plate,  so  hf+  ■  c^  +e,  where  h^+  is  the 
±  thrust  plate  gap. 

Let  f /M  *  a  ,  the  external  acceleration,  then  (24)  and  (25)  give 
e 

2  2 

-md  e/dt  **  mag  -  F  (26) 

Let  e  “  ®/c0>  in  nondimens icnal  form  (26)  is  . 

d2e/dt2  =  (-ag  +  F/m)/c  v2,  (27) 

.-nd  it  is  further  convenient  to  put  (27)  in  system  form 
de^/dT  -  e2 


de,/dT  »  (-a  +  F/m)/c  v' 


(28) 


16. 


where  s,  =  e/c  and  e„  ■  (de/dt)/c  v. 

I  o  <s  o 

Different  methods  are  used  to  integrate  (28)  depending  on  the  truncation 
error  of  the  scheme  used  to  integrate  the  Reynolds  equation. 
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The  spatial  difference  approximations  to  (20)  are  obtained  by  conventional 
divergence  (or  integral)  methods.  An  arbitrary  mesh  is  used  with  discontinuities 
in  any  bearing  parameter  only  being  allowed  to  occur  at  a  mesh  point.  Mesh  points 
in  the  direction  of  increasing  S  are  enumerated  by  the  index  J.  Define 
AEj  =  Sj  +  ^  -S  ,  and  define  the  fundamental  interval  i(J)  to  be 

KJ)  =  [Sj  -  .5  ASj  r  Sj  +  .5  ASj]  (29) 

which  has  length  .5(ASJ_1  +  ASj).  Integrate  ( 2  Q)  over  1(J)  to  obtain 


r  4.  r-L  T  -Yc  T  e  *  A-~ -  R  (k  H  ||  +  P 

S,j  t  *5  S,J+  S,J"  S, J  -.5  j  \  0)  r  BT  BT  / 


Since  a  L  _  except  at  the  chamfer,  by  (18)  (30)  becomes 

b  j  JT  b  j  J" 

Y-  +  s  -Y<?  s  a  A  ^-7 -  r(k  hJ|  +  P  )ds 

S,J  +  .5  S, J  -.5  j  cu ^  +iu2  \  u)  raT  iT  / 

_Al°  _2v _  R  =  ,p  8ik  +  £ 

Aro  U) u) 2  Rj  h  (P  8T  Ch  &T) 
where  the  second  term  on  the  right  appears  only  at  thechamfer. 


De  f ine 


a  -  [PhJKjR] 


b  =  [AcK.R  cos0  P] 

0  4 

g  -A  — ^ -  [R  K  H  3 

oj  ^  +  oj  2  »  r 

c  =  A  — Q —  [rp  &h  /St] 
W1  +  w2  r 

V  ’  <pj  +  i  -pj>/1sj 


The  basic  approximation  to  be  used  to  (31)  is 
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3J  +  .5AJP  'bJ  +  .5  ~aj  -  ,5AJ  -  1P  +  bJ  -  .5 


=  -5  (ASJ  gJ+  +  ASJ-1  8J-) 


Pn+1  .  pn 

\J _ \j 

A  T 


.5  («Sj. 


CJ+  +  “j-l'j- 
P"+1  -  Prf 


+  A  ^^7  RJ  £h  (PJ  t8Hr,J+/ai  +  *  S  -JTf  ~) 


where  superscript  n  refers  to  the  n-th  time  step.  Subscript  J  refers  to  the 
J-th  spatial  mesh  point.  All  quantities  defined  at  the  spatial  half  step 
that  cannot  be  analytically  directly  determined  are  obtained  by  linear 
interpolation,  e.g.  PJ+  ^  =  -5(Pj  +  Pj+p.  Quantities  such  as  gj+  and  g^_ 
are  respectively  the  right  and  left  limits  at  the  J-th  mesh  point  of  g.  The 
replacement  of  dH  /§T  by  dH  /dT  +  dH  /dT  is  valid  since  one  term  or  the 
other  is  0  at  the  chamfer. 


(33) 
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APPENDIX  4.  Geometry  of  Spool  Bearing 


The  pertinent  parameters  for  the  spool  bearing  examples  given 
are  summarized  in  Table  4. 


TABLE  4.  SPOOL  BEARING  PARAMETERS 


Parameter 

Thrust 

Journal 

a 

.5 

.5 

0  (degrees) 

166 

30 

6  (in.) 

1 36 ( 10) -6 

136(10) 

hr  (in.) 

90(10)'6 

90(10) 

Y 

.65 

.40 

ro  (in<) 

.52 

- 

r1  (in.) 

.23 

- 

L  (in.) 

- 

1.8 

r.  (in.) 

- 

.23 

u2  (RPM) 

0 

24,000 

U1  (RPM) 

24,000 

0 

Other  quantities  are  p. 

d 

and  m  *  49  gm. 


«  29.4  lb. /in.2. 


u  =  3.02  (10)"9  lb.  sec/in.2 


NOMENCLATURE 


Upper  Case  English 


F 

H 

H, 

H 

J 

J± 


r± 


K 

K, 


L 

M 

P 

R 


cf/co*  ^dimensional  chamfer  height 

gas  film  force  in  axial  direction  (lb.) 

nondimensaional  bearing  gap  for  plain  journal  bearing 

hr/cQ,  nondimensional  bearing  gap 

at  positive  or  negative  thrust  plate 

subscript  indicating  J-th  spatial  mesh  point 

subscript  indicating  right  or  left  limit  at  J-th 
mesh  point 

?  ?  1/2 

(U  +  V  )  ,  total  bearing  stiffness  (lb. /in.) 

Whipple  coefficient  [See  (15)] 

H 

II 

length  of  journal  bearing  (in.) 
mass  of  stator  (lb.  sec.  /in.) 
p/p  .  nondimensional  pressure 

a 

r/r  ,  nondimensional  radial  distance  from  axis 
of  symmetry. 


N-2 


S 

T 

T' 

U 

V 


s/rQ,  nondimensional  meridianal  coordinate 

vt,  nondimensional  time 

T/(2ir/288),  abscissa  of  figures 

in-phase  bearing  stiffness 

out-of-phase  bearing  stiffness 

constant  term  of  V  =  VQ  +  V^/N* 

coefficient  of  linear  term  of  V  =  V  +  V,/N* 

o  l 

ratio  of  grooved  region  to  total  bearing  extent 
in  the  meridianal  direction 


coefficient  of  sin(T)  term  [See  (5)] 


N 


coefficient  for  difference  approximation  [See  (32)] 

nominal  clearance  used  to  nondimensional ize  quantities 
(in.) 

chamfer  height  (in.) 

external  force  on  bearing  (lb.) 

coefficient  for  difference  approximation  [See  (32)] 

bearing  clearance  (in.) 

hr  at  positive  or  negative  thrust  plate 

mesh  interval  from  midpoint  to  midpoint  [See  (29)] 

rotor  mass  (lb.  sec.  /in.) 

superscript  designating  time  [e.g.,  Pn  =  P(haT)] 
o 

pressure  (lb. /in.  ) 

ambient  pressure  (lb. /in.  ) 

radial  distance  from  axis  of  symmetry  (in.). 

outer  radius  of  thrust  plate  (in.) 

inner  radius  of  thrust  plate  (in.) 

journal  radius  (in.) 

meridianal  coordinate  (in.) 


time  (sec.) 


displacement  of  stator  (in.) 


=  displacement  of  rotor  (in.) 

Upper  Case  Greek 

r  =  (hr  +  <5)/hr,  groove  depth  ratio 

A  =  forward  difference  with  respect  to  S,  e.g., 

J  4jp  ■  <PJ+1  -  pj)/as0 

aSj  =  mesh  spacing  at  J-th  mesh  point,  AS^  =  SJ+1  -  Sj 

aT  =  nondimensional  time  step 

A  =  compressibility  number  [See  (15)] 

A.  =  Whipple  coefficient  [See  (15)] 

0 

¥<.  =  nondimensional  mass  flux  [See  (14)] 

¥  =  nondimensional  mass  flux  flowing  into  chamfer 

S|in 

•  .  =  nondimensional  mass  flux  flowing  out  of  chamfer 

S|out 

Lower  Class  Greek 

a  =  ratio  of  groove  width  to  sum  of  groove  and  ridge 

width 

8  =  spiral  groove  angle 


6 


groove  depth  (in.) 


N-5 


dimensional  axial  eccentricity  (in.)  [See  (25)] 
e/cQ ,  nondimensional  eccentricity 
e/c0,  nondimensional  eccentricity  (Same  as  e) 
(d€/dt)/cQv,  nondimensional  eccentricity  velocity 

p 

viscosity  (lb. sec. /in.  ) 

frequency  (radians/sec.) 

arctan  (-V/U),  phase  angle  (radians) 

smooth  surface  rotational  speed  (radians/sec.) 

grooved  surface  rotational  speed  (radians/sec.) 


subscript  used  to  denote  positive  and  negative  thrust 
plate  or  when  used  as  J±  indicates  right  or  left  limit 
at  J-th  mesh  point 


ECCENTRICITY 
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Figure  2 


Implicit  Scheme  Axial  Dynamics  for  Very  Fine 


ECCENTRICITY 


Figure  3  -  Extrapolated  Crank-Nicolson  Scheme  at  Coarse 

Time  Step 
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